home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 40 / Amiga Format CD40 (1999-05-11)(Future Publishing)(GB)(Track 1 of 3)[!][issue 1999-06].iso / -readerstuff- / paul_qureshi / source / tritri.c < prev    next >
C/C++ Source or Header  |  1999-03-27  |  10KB  |  292 lines
  1. /* Triangle/triangle intersection test routine,
  2.  * by Tomas Moller, 1997.
  3.  * See article "A Fast Triangle-Triangle Intersection Test",
  4.  * Journal of Graphics Tools, 2(2), 1997
  5.  *
  6.  * int tri_tri_intersect(float V0[3],float V1[3],float V2[3],
  7.  *                         float U0[3],float U1[3],float U2[3])
  8.  *
  9.  * parameters: vertices of triangle 1: V0,V1,V2
  10.  *             vertices of triangle 2: U0,U1,U2
  11.  * result    : returns 1 if the triangles intersect, otherwise 0
  12.  *
  13.  */
  14.  
  15. #include <math.h>
  16.  
  17.  
  18. /* if USE_EPSILON_TEST is true then we do a check: 
  19.          if |dv|<EPSILON then dv=0.0;
  20.    else no check is done (which is less robust)
  21. */
  22. #define USE_EPSILON_TEST TRUE
  23. #define EPSILON 0.000001
  24.  
  25.  
  26. /* some macros */
  27. #define CROSS(dest,v1,v2)                      \
  28.               dest[0]=v1[1]*v2[2]-v1[2]*v2[1]; \
  29.               dest[1]=v1[2]*v2[0]-v1[0]*v2[2]; \
  30.               dest[2]=v1[0]*v2[1]-v1[1]*v2[0];
  31.  
  32. #define DOT(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2])
  33.  
  34. #define SUB(dest,v1,v2)          \
  35.             dest[0]=v1[0]-v2[0]; \
  36.             dest[1]=v1[1]-v2[1]; \
  37.             dest[2]=v1[2]-v2[2]; 
  38.  
  39. /* sort so that a<=b */
  40. #define SORT(a,b)       \
  41.              if(a>b)    \
  42.              {          \
  43.                float c; \
  44.                c=a;     \
  45.                a=b;     \
  46.                b=c;     \
  47.              }
  48.  
  49. #define ISECT(VV0,VV1,VV2,D0,D1,D2,isect0,isect1) \
  50.               isect0=VV0+(VV1-VV0)*D0/(D0-D1);    \
  51.               isect1=VV0+(VV2-VV0)*D0/(D0-D2);
  52.  
  53.  
  54. #define COMPUTE_INTERVALS(VV0,VV1,VV2,D0,D1,D2,D0D1,D0D2,isect0,isect1) \
  55.   if(D0D1>0.0f)                                         \
  56.   {                                                     \
  57.     /* here we know that D0D2<=0.0 */                   \
  58.     /* that is D0, D1 are on the same side, D2 on the other or on the plane */ \
  59.     ISECT(VV2,VV0,VV1,D2,D0,D1,isect0,isect1);          \
  60.   }                                                     \
  61.   else if(D0D2>0.0f)                                    \
  62.   {                                                     \
  63.     /* here we know that d0d1<=0.0 */                   \
  64.     ISECT(VV1,VV0,VV2,D1,D0,D2,isect0,isect1);          \
  65.   }                                                     \
  66.   else if(D1*D2>0.0f || D0!=0.0f)                       \
  67.   {                                                     \
  68.     /* here we know that d0d1<=0.0 or that D0!=0.0 */   \
  69.     ISECT(VV0,VV1,VV2,D0,D1,D2,isect0,isect1);          \
  70.   }                                                     \
  71.   else if(D1!=0.0f)                                     \
  72.   {                                                     \
  73.     ISECT(VV1,VV0,VV2,D1,D0,D2,isect0,isect1);          \
  74.   }                                                     \
  75.   else if(D2!=0.0f)                                     \
  76.   {                                                     \
  77.     ISECT(VV2,VV0,VV1,D2,D0,D1,isect0,isect1);          \
  78.   }                                                     \
  79.   else                                                  \
  80.   {                                                     \
  81.     /* triangles are coplanar */                        \
  82.     return coplanar_tri_tri(N1,V0,V1,V2,U0,U1,U2);      \
  83.   }
  84.  
  85.  
  86.  
  87. /* this edge to edge test is based on Franlin Antonio's gem:
  88.    "Faster Line Segment Intersection", in Graphics Gems III,
  89.    pp. 199-202 */ 
  90. #define EDGE_EDGE_TEST(V0,U0,U1)                      \
  91.   Bx=U0[i0]-U1[i0];                                   \
  92.   By=U0[i1]-U1[i1];                                   \
  93.   Cx=V0[i0]-U0[i0];                                   \
  94.   Cy=V0[i1]-U0[i1];                                   \
  95.   f=Ay*Bx-Ax*By;                                      \
  96.   d=By*Cx-Bx*Cy;                                      \
  97.   if((f>0 && d>=0 && d<=f) || (f<0 && d<=0 && d>=f))  \
  98.   {                                                   \
  99.     e=Ax*Cy-Ay*Cx;                                    \
  100.     if(f>0)                                           \
  101.     {                                                 \
  102.       if(e>=0 && e<=f) return 1;                      \
  103.     }                                                 \
  104.     else                                              \
  105.     {                                                 \
  106.       if(e<=0 && e>=f) return 1;                      \
  107.     }                                                 \
  108.   }                                
  109.  
  110. #define EDGE_AGAINST_TRI_EDGES(V0,V1,U0,U1,U2) \
  111. {                                              \
  112.   float Ax,Ay,Bx,By,Cx,Cy,e,d,f;               \
  113.   Ax=V1[i0]-V0[i0];                            \
  114.   Ay=V1[i1]-V0[i1];                            \
  115.   /* test edge U0,U1 against V0,V1 */          \
  116.   EDGE_EDGE_TEST(V0,U0,U1);                    \
  117.   /* test edge U1,U2 against V0,V1 */          \
  118.   EDGE_EDGE_TEST(V0,U1,U2);                    \
  119.   /* test edge U2,U1 against V0,V1 */          \
  120.   EDGE_EDGE_TEST(V0,U2,U0);                    \
  121. }
  122.  
  123. #define POINT_IN_TRI(V0,U0,U1,U2)           \
  124. {                                           \
  125.   float a,b,c,d0,d1,d2;                     \
  126.   /* is T1 completly inside T2? */          \
  127.   /* check if V0 is inside tri(U0,U1,U2) */ \
  128.   a=U1[i1]-U0[i1];                          \
  129.   b=-(U1[i0]-U0[i0]);                       \
  130.   c=-a*U0[i0]-b*U0[i1];                     \
  131.   d0=a*V0[i0]+b*V0[i1]+c;                   \
  132.                                             \
  133.   a=U2[i1]-U1[i1];                          \
  134.   b=-(U2[i0]-U1[i0]);                       \
  135.   c=-a*U1[i0]-b*U1[i1];                     \
  136.   d1=a*V0[i0]+b*V0[i1]+c;                   \
  137.                                             \
  138.   a=U0[i1]-U2[i1];                          \
  139.   b=-(U0[i0]-U2[i0]);                       \
  140.   c=-a*U2[i0]-b*U2[i1];                     \
  141.   d2=a*V0[i0]+b*V0[i1]+c;                   \
  142.   if(d0*d1>0.0)                             \
  143.   {                                         \
  144.     if(d0*d2>0.0) return 1;                 \
  145.   }                                         \
  146. }
  147.  
  148. int coplanar_tri_tri(float N[3],float V0[3],float V1[3],float V2[3],
  149.                      float U0[3],float U1[3],float U2[3])
  150. {
  151.    float A[3];
  152.    short i0,i1;
  153.    /* first project onto an axis-aligned plane, that maximizes the area */
  154.    /* of the triangles, compute indices: i0,i1. */
  155.    A[0]=fabs(N[0]);
  156.    A[1]=fabs(N[1]);
  157.    A[2]=fabs(N[2]);
  158.    if(A[0]>A[1])
  159.    {
  160.       if(A[0]>A[2])
  161.       {
  162.           i0=1;      /* A[0] is greatest */
  163.           i1=2;
  164.       }
  165.       else
  166.       {
  167.           i0=0;      /* A[2] is greatest */
  168.           i1=1;
  169.       }
  170.    }
  171.    else   /* A[0]<=A[1] */
  172.    {
  173.       if(A[2]>A[1])
  174.       {
  175.           i0=0;      /* A[2] is greatest */
  176.           i1=1;
  177.       }
  178.       else
  179.       {
  180.           i0=0;      /* A[1] is greatest */
  181.           i1=2;
  182.       }
  183.     }               
  184.                 
  185.     /* test all edges of triangle 1 against the edges of triangle 2 */
  186.     EDGE_AGAINST_TRI_EDGES(V0,V1,U0,U1,U2);
  187.     EDGE_AGAINST_TRI_EDGES(V1,V2,U0,U1,U2);
  188.     EDGE_AGAINST_TRI_EDGES(V2,V0,U0,U1,U2);
  189.                 
  190.     /* finally, test if tri1 is totally contained in tri2 or vice versa */
  191.     POINT_IN_TRI(V0,U0,U1,U2);
  192.     POINT_IN_TRI(U0,V0,V1,V2);
  193.  
  194.     return 0;
  195. }
  196.  
  197.  
  198. int tri_tri_intersect(float V0[3],float V1[3],float V2[3],
  199.                       float U0[3],float U1[3],float U2[3])
  200. {
  201.   float E1[3],E2[3];
  202.   float N1[3],N2[3],d1,d2;
  203.   float du0,du1,du2,dv0,dv1,dv2;
  204.   float D[3];
  205.   float isect1[2], isect2[2];
  206.   float du0du1,du0du2,dv0dv1,dv0dv2;
  207.   short index;
  208.   float vp0,vp1,vp2;
  209.   float up0,up1,up2;
  210.   float b,c,max;
  211.  
  212.   /* compute plane equation of triangle(V0,V1,V2) */
  213.   SUB(E1,V1,V0);
  214.   SUB(E2,V2,V0);
  215.   CROSS(N1,E1,E2);
  216.   d1=-DOT(N1,V0);
  217.   /* plane equation 1: N1.X+d1=0 */
  218.  
  219.   /* put U0,U1,U2 into plane equation 1 to compute signed distances to the plane*/
  220.   du0=DOT(N1,U0)+d1;
  221.   du1=DOT(N1,U1)+d1;
  222.   du2=DOT(N1,U2)+d1;
  223.  
  224.   /* coplanarity robustness check */
  225. #if USE_EPSILON_TEST==TRUE
  226.   if(fabs(du0)<EPSILON) du0=0.0;
  227.   if(fabs(du1)<EPSILON) du1=0.0;
  228.   if(fabs(du2)<EPSILON) du2=0.0;
  229. #endif
  230.   du0du1=du0*du1;
  231.   du0du2=du0*du2;
  232.  
  233.   if(du0du1>0.0f && du0du2>0.0f) /* same sign on all of them + not equal 0 ? */
  234.     return 0;                    /* no intersection occurs */
  235.  
  236.   /* compute plane of triangle (U0,U1,U2) */
  237.   SUB(E1,U1,U0);
  238.   SUB(E2,U2,U0);
  239.   CROSS(N2,E1,E2);
  240.   d2=-DOT(N2,U0);
  241.   /* plane equation 2: N2.X+d2=0 */
  242.  
  243.   /* put V0,V1,V2 into plane equation 2 */
  244.   dv0=DOT(N2,V0)+d2;
  245.   dv1=DOT(N2,V1)+d2;
  246.   dv2=DOT(N2,V2)+d2;
  247.  
  248. #if USE_EPSILON_TEST==TRUE
  249.   if(fabs(dv0)<EPSILON) dv0=0.0;
  250.   if(fabs(dv1)<EPSILON) dv1=0.0;
  251.   if(fabs(dv2)<EPSILON) dv2=0.0;
  252. #endif
  253.  
  254.   dv0dv1=dv0*dv1;
  255.   dv0dv2=dv0*dv2;
  256.         
  257.   if(dv0dv1>0.0f && dv0dv2>0.0f) /* same sign on all of them + not equal 0 ? */
  258.     return 0;                    /* no intersection occurs */
  259.  
  260.   /* compute direction of intersection line */
  261.   CROSS(D,N1,N2);
  262.  
  263.   /* compute and index to the largest component of D */
  264.   max=fabs(D[0]);
  265.   index=0;
  266.   b=fabs(D[1]);
  267.   c=fabs(D[2]);
  268.   if(b>max) max=b,index=1;
  269.   if(c>max) max=c,index=2;
  270.  
  271.         /* this is the simplified projection onto L*/
  272.         vp0=V0[index];
  273.         vp1=V1[index];
  274.         vp2=V2[index];
  275.  
  276.         up0=U0[index];
  277.         up1=U1[index];
  278.         up2=U2[index];
  279.  
  280.  /* compute interval for triangle 1 */
  281.   COMPUTE_INTERVALS(vp0,vp1,vp2,dv0,dv1,dv2,dv0dv1,dv0dv2,isect1[0],isect1[1]);
  282.  
  283.  /* compute interval for triangle 2 */
  284.   COMPUTE_INTERVALS(up0,up1,up2,du0,du1,du2,du0du1,du0du2,isect2[0],isect2[1]);
  285.  
  286.   SORT(isect1[0],isect1[1]);
  287.   SORT(isect2[0],isect2[1]);
  288.  
  289.   if(isect1[1]<isect2[0] || isect2[1]<isect1[0]) return 0;
  290.   return 1;
  291. }
  292.